function [C_th] = theo(rho,var_e,var_eta,var_z0,var_a) 

T=29;

% this refers to whether one thinks of first period as 0 or 1 (i.e. with this,
% process "starts" one period later)
% results are very similiar for both cases
var_z0 = var_eta;

 % theoretical C matrix
 C_th = NaN(T,T);
 
 %% diagonal entries
 C_th(1,1) = rho^2*var_z0+ var_eta+var_e+var_a;
 
 C_th(2,2) = rho^(2*2)*var_z0+rho^(2*1)*var_eta +var_eta+ var_e+var_a;
 
 C_th(3,3) = rho^(2*3)*var_z0+rho^(2*2)*var_eta+rho^(2*1)*var_eta+var_eta + var_e+var_a;
 
 C_th(4,4) = rho^(2*4)*var_z0+rho^(2*3)*var_eta+rho^(2*2)*var_eta+rho^(2*1)*var_eta+var_eta + var_e+var_a;
 
 C_th(5,5) = rho^(2*5)*var_z0 +rho^(2*4)*var_eta+rho^(2*3)*var_eta+...
             rho^(2*2)*var_eta+rho^(2*1)*var_eta+var_eta + var_e+var_a;
         
 C_th(6,6) = rho^(2*6)*var_z0+rho^(2*5)*var_eta+rho^(2*4)*var_eta+rho^(2*3)*var_eta+rho^(2*2)*var_eta+...
             rho^(2*1)*var_eta+var_eta + var_e+var_a;
         
 C_th(7,7) = rho^(2*7)*var_z0+rho^(2*6)*var_eta+rho^(2*5)*var_eta+rho^(2*4)*var_eta+rho^(2*3)*var_eta+...
             rho^(2*2)*var_eta+rho^(2*1)*var_eta+var_eta + var_e+var_a;
         
 C_th(8,8) = rho^(2*8)*var_z0+rho^(2*7)*var_eta+rho^(2*6)*var_eta+rho^(2*5)*var_eta+rho^(2*4)*var_eta+...
             rho^(2*3)*var_eta+rho^(2*2)*var_eta+rho^(2*1)*var_eta+var_eta + var_e+var_a;
 
 C_th(9,9) = rho^(2*9)*var_z0+rho^(2*8)*var_eta+rho^(2*7)*var_eta+rho^(2*6)*var_eta+rho^(2*5)*var_eta+...
             rho^(2*4)*var_eta+rho^(2*3)*var_eta+rho^(2*2)*var_eta+rho^(2*1)*var_eta+var_eta + var_e+var_a;
 
 C_th(10,10) = rho^(2*10)*var_z0+rho^(2*9)*var_eta+rho^(2*8)*var_eta+rho^(2*7)*var_eta+rho^(2*6)*var_eta+...
               rho^(2*5)*var_eta+rho^(2*4)*var_eta+rho^(2*3)*var_eta+rho^(2*2)*var_eta+rho^(2*1)*var_eta+...
               var_eta + var_e+var_a;

 C_th(11,11) = rho^(2*11)*var_z0+rho^(2*10)*var_eta+rho^(2*9)*var_eta+rho^(2*8)*var_eta+rho^(2*7)*var_eta+...
               rho^(2*6)*var_eta+rho^(2*5)*var_eta+rho^(2*4)*var_eta+rho^(2*3)*var_eta+rho^(2*2)*var_eta+...
               rho^(2*1)*var_eta+var_eta + var_e+var_a;

 C_th(12,12) = rho^(2*12)*var_z0+rho^(2*11)*var_eta+rho^(2*10)*var_eta+rho^(2*9)*var_eta+rho^(2*8)*var_eta+...
               rho^(2*7) *var_eta+rho^(2*6)*var_eta+rho^(2*5)*var_eta+rho^(2*4)*var_eta+rho^(2*3)*var_eta+...
               rho^(2*2) *var_eta+rho^(2*1)*var_eta+var_eta + var_e+var_a;
 
 C_th(13,13) = rho^(2*13)*var_z0 +rho^(2*12)*var_eta+rho^(2*11)*var_eta+rho^(2*10)*var_eta+...
               rho^(2*9) *var_eta+rho^(2*8)*var_eta+rho^(2*7)*var_eta+rho^(2*6)*var_eta+...
               rho^(2*5)*var_eta +rho^(2*4)*var_eta+rho^(2*3)*var_eta+rho^(2*2)*var_eta+...
               rho^(2*1)*var_eta +var_eta + var_e+var_a;

 C_th(14,14) = rho^(2*14)*var_z0+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+rho^(2*10)*var_eta+...
               rho^(2*9) *var_eta+rho^(2*8)*var_eta+rho^(2*7) *var_eta+rho^(2*6)*var_eta+rho^(2*5)*var_eta+...
               rho^(2*4) *var_eta+rho^(2*3)*var_eta+rho^(2*2) *var_eta+rho^(2*1)*var_eta+var_eta + var_e+var_a;

 C_th(15,15) = rho^(2*15)*var_z0+rho^(2*14) *var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+...
               rho^(2*11)*var_eta+rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8)*var_eta+...
               rho^(2*7) *var_eta+rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4)*var_eta+...
               rho^(2*3) *var_eta+rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a;

 C_th(16,16) = rho^(2*16)*var_z0 +rho^(2*15)*var_eta+rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+...
               rho^(2*11)*var_eta+rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8)*var_eta+...
               rho^(2*7) *var_eta+rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4)*var_eta+...
               rho^(2*3) *var_eta+rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a;

 C_th(17,17) = rho^(2*17)*var_z0 +rho^(2*16)*var_eta+rho^(2*15)*var_eta+rho^(2*14)*var_eta+...
               rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+rho^(2*10)*var_eta+...
               rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+rho^(2*6) *var_eta+...
               rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+rho^(2*2) *var_eta+...
               rho^(2*1) *var_eta+var_eta + var_e+var_a; 

 C_th(18,18) = rho^(2*18)*var_z0 +rho^(2*17)*var_eta+rho^(2*16)*var_eta+rho^(2*15)*var_eta+...
               rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+...
               rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+...
               rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+...
               rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a;

 C_th(19,19) = rho^(2*19)*var_z0+rho^(2*18) *var_eta+rho^(2*17)*var_eta+rho^(2*16)*var_eta+rho^(2*15)*var_eta+...
               rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+...
               rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+...
               rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+rho^(2*2)*var_eta+...
               rho^(2*1) *var_eta+var_eta + var_e+var_a;
 
 C_th(20,20) = rho^(2*20)*var_z0 +rho^(2*19)*var_eta+rho^(2*18)*var_eta+rho^(2*17)*var_eta+rho^(2*16)*var_eta+...
               rho^(2*15)*var_eta+rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+...
               rho^(2*11)*var_eta+rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8)*var_eta+...
               rho^(2*7) *var_eta+rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4)*var_eta+...
               rho^(2*3) *var_eta+rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a;
 
 C_th(21,21) = rho^(2*21)*var_z0 +rho^(2*20)*var_eta+rho^(2*19)*var_eta+rho^(2*18)*var_eta+rho^(2*17)*var_eta+rho^(2*16)*var_eta+...
               rho^(2*15)*var_eta+rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+...
               rho^(2*11)*var_eta+rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+...
               rho^(2*7) *var_eta+rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+...
               rho^(2*3)*var_eta+rho^(2*2)*var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a;
 
 C_th(22,22) = rho^(2*22)*var_z0 +rho^(2*21)*var_eta+rho^(2*20)*var_eta+rho^(2*19)*var_eta+...
               rho^(2*18)*var_eta+rho^(2*17)*var_eta+rho^(2*16)*var_eta+rho^(2*15)*var_eta+...
               rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+...
               rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+...
               rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+...
               rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a ;
           
 C_th(23,23) = rho^(2*23)*var_z0 +rho^(2*22)*var_eta+rho^(2*21)*var_eta+rho^(2*20)*var_eta+rho^(2*19)*var_eta+...
               rho^(2*18)*var_eta+rho^(2*17)*var_eta+rho^(2*16)*var_eta+rho^(2*15)*var_eta+...
               rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+...
               rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+...
               rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+...
               rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a ;
           
 C_th(24,24) = rho^(2*24)*var_z0 +rho^(2*23)*var_eta+rho^(2*22)*var_eta+rho^(2*21)*var_eta+...
               rho^(2*20)*var_eta+rho^(2*19)*var_eta+rho^(2*18)*var_eta+rho^(2*17)*var_eta+...
               rho^(2*16)*var_eta+rho^(2*15)*var_eta+...
               rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+...
               rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+...
               rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+...
               rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a ;

 C_th(25,25) = rho^(2*25)*var_z0 +rho^(2*24)*var_eta +rho^(2*23)*var_eta+rho^(2*22)*var_eta+rho^(2*21)*var_eta+...
               rho^(2*20)*var_eta+rho^(2*19)*var_eta+rho^(2*18)*var_eta+rho^(2*17)*var_eta+...
               rho^(2*16)*var_eta+rho^(2*15)*var_eta+...
               rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+...
               rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+...
               rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+...
               rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a ;
           
 C_th(26,26) = rho^(2*26)*var_z0 +rho^(2*25)*var_eta +rho^(2*24)*var_eta +rho^(2*23)*var_eta+...
               rho^(2*22)*var_eta+rho^(2*21)*var_eta+...
               rho^(2*20)*var_eta+rho^(2*19)*var_eta+rho^(2*18)*var_eta+rho^(2*17)*var_eta+...
               rho^(2*16)*var_eta+rho^(2*15)*var_eta+...
               rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+...
               rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+...
               rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+...
               rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a ;
           
 C_th(27,27) = rho^(2*27)*var_z0 +rho^(2*26)*var_eta +rho^(2*25)*var_eta +rho^(2*24)*var_eta+...
               rho^(2*23)*var_eta+rho^(2*22)*var_eta+rho^(2*21)*var_eta+...
               rho^(2*20)*var_eta+rho^(2*19)*var_eta+rho^(2*18)*var_eta+rho^(2*17)*var_eta+...
               rho^(2*16)*var_eta+rho^(2*15)*var_eta+...
               rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+...
               rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+...
               rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+...
               rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a ;
         
 C_th(28,28) = rho^(2*28)*var_z0 +rho^(2*27)*var_eta +rho^(2*26)*var_eta +rho^(2*25)*var_eta+...
               rho^(2*24)*var_eta+rho^(2*23)*var_eta+rho^(2*22)*var_eta+rho^(2*21)*var_eta+...
               rho^(2*20)*var_eta+rho^(2*19)*var_eta+rho^(2*18)*var_eta+rho^(2*17)*var_eta+...
               rho^(2*16)*var_eta+rho^(2*15)*var_eta+...
               rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+...
               rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+...
               rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+...
               rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a ;
           
 C_th(29,29) = rho^(2*29)*var_z0 +rho^(2*28)*var_eta +rho^(2*27)*var_eta +rho^(2*26)*var_eta+...
               rho^(2*25)*var_eta+rho^(2*24)*var_eta+rho^(2*23)*var_eta+rho^(2*22)*var_eta+rho^(2*21)*var_eta+...
               rho^(2*20)*var_eta+rho^(2*19)*var_eta+rho^(2*18)*var_eta+rho^(2*17)*var_eta+...
               rho^(2*16)*var_eta+rho^(2*15)*var_eta+...
               rho^(2*14)*var_eta+rho^(2*13)*var_eta+rho^(2*12)*var_eta+rho^(2*11)*var_eta+...
               rho^(2*10)*var_eta+rho^(2*9) *var_eta+rho^(2*8) *var_eta+rho^(2*7) *var_eta+...
               rho^(2*6) *var_eta+rho^(2*5) *var_eta+rho^(2*4) *var_eta+rho^(2*3) *var_eta+...
               rho^(2*2) *var_eta+rho^(2*1) *var_eta+var_eta + var_e+var_a ;
           

%% offdiagnoal entries 
 for i = 1:T
    for j = i+1:T
 C_th(i,j) = var_a+rho^(j-i)*(C_th(i,i)-var_e);
    end
 end
 
 %% stack C matrix

C_th=reshape(C_th,[T*T,1]);  
      id1 = find(isnan(C_th));
        C_th([id1]) = [];
end